# draw program for human


data <- read.table(file.choose(), header=F)


#===============================================================
# change angle
#===============================================================
Transform <- function(mat){
	d = dim(mat)[1]
	n = dim(mat)[2]
	mat.rev = t(mat[d+1-c(1:d), ])
	mat.rev
}

c <- colorRampPalette(c("#00007F", "blue", "#007FFF", 
	"cyan","#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"), bias=1.5)( 100 )
c <- colorRampPalette(c("#00007F", "blue", "#007FFF", 
	"cyan","#7FFF7F", "yellow", "#FF7F00", "red", "#7F0000"))( 100 )
#c <- colorRampPalette(c("#FFFFD4", "#FED98E", "#FE9929", "#D95F0E", "#993404"),space = "Lab")(100)	# yellow to orrange


#===============================================================
# Extract part of the graph
#===============================================================

mat <- matrix(data$V3, nrow=sqrt(length(data$V3)), byrow=T)

index <- intersect(grep("chr14:",data$V1), grep("chr14:", data$V2))
mat <- matrix(data$V3[index], nrow=sqrt(length(data$V3[index])), byrow=T)
mat.rev <- Transform(mat)

#===============================================================
# Modify value
#===============================================================
### setting maximum value
#maxVal <- max(mat)
maxVal <- 5
mat.conv <- ifelse(mat.rev > maxVal, maxVal, mat.rev)
minVal <- 0.2
mat.conv <- ifelse(mat.conv < minVal, minVal, mat.conv)


#===============================================================
# plot start
#---------------------------------------------------------------
title <- 'HiC human contact map (correlation)'
subTitle <- '1Mb resolution(chr14)'

li <- layout(
	matrix(
		c(rep(1,10), 2
		),
		byrow=T, ncol=1, nrow=11
	)
)


par(oma=c(3,4,6,4), mar=c(0,0,0,0))

# image 1 :association map
image(mat.conv, col=c, axes=F)

# image 2 :color chart
library(fields)
plot(c(0,1),c(0,1),type='n', xlab='', ylab='', main="", axes=F)
image.plot(mat.conv, col=c, horizontal=T, axes=F, legend.only=T, legend.width=5)

mtext(side = 3, outer=T, text=title, cex=1.2, line=3)
mtext(side = 3, outer=T, text=subTitle, cex=0.8, line=1.5)

#---------------------------------------------------------------
# end
#===============================================================







#===============================================================
# plot with FISH locations 40kb resolution map
#---------------------------------------------------------------

title <- 'After 5 SCN normalization 1%'
subTitle <- 'remove duplicate, max=0.15, normalization of raw score'

li <- layout(
	matrix(
		c(
			rep(1,13), 5,
			rep(c(rep(3,13),2),13),
			rep(c(rep(4,13),6),2)
		),
		byrow=T, ncol=14, nrow=16
	)
)


par(oma=c(3,6,6,4), mar=c(0,0,0,0))

# image 1 :chromosome
plot(image.h)

# image 2 :chromosome
plot(image.v)

# image 3 :association map
image(mat.conv, col=c, axes=F)


FISH <- function(x,y,name, upper=12, adj=0.5){
	par(new=T, xaxs='i', yaxs='i')
	plot(x,316+1-y, cex=2, xlim=c(1,316), ylim=c(1,316), axes=F, col='white', lwd=3)
	text(x,316+1-y+upper, cex=1.3, name, col='white', adj=adj, font=2)
}
FISH(220,218,'c609-c1778')
FISH(20,18,'c4G8-c22F3', upper=-8)
FISH(294,290,'c162-c622')
FISH(182,38,'c343-c5D6')
FISH(45,38,'c17A5-c5D6', upper=-8)
FISH(195,187,'c17G9-c3H7')
FISH(122,110,'c20G4-c23D3')
FISH(22,5,'c23G3-c13G6', upper=1.5,adj=0)	#
FISH(38,20,'c5D6-c4G8', upper=1.5,adj=0)
FISH(76,38,'c1786-c5D6', upper=1.5, adj=0)
FISH(301,256,'c1223-c757', upper=1.5, adj=1)
FISH(301,230,'c1223-c16D10', upper=1.5, adj=1)
FISH(302,21,'c737-c2G11',upper=1.5, adj=1)
FISH(211,25,'c821-c19C7')
FISH(195,107,'c17G9-c2C4')
FISH(302,69,'c6F6-c737', upper=1.5, adj=1)
FISH(302,49,'c737-c3H1', upper=1.5, adj=1)
FISH(290,218,'c1778-c622')


# image 4 :color chart
library(fields)
plot(c(0,1),c(0,1),type='n', xlab='', ylab='', main="", axes=F)
image.plot(mat.conv, col=c, horizontal=T, axes=F, legend.only=T, legend.width=5)

# image 5-6:dummy
plot(c(0,1),c(0,1),type='n', xlab='', ylab='', main="", axes=F)
plot(c(0,1),c(0,1),type='n', xlab='', ylab='', main="", axes=F)

mtext(side = 3, outer=T, text=title, cex=1.2, line=3)
mtext(side = 3, outer=T, text=subTitle, cex=0.8, line=1.5)

#---------------------------------------------------------------
# end
#===============================================================


